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We study the properties of a nano-electromechanical system in the coherent regime, where the 
electronic and vibrational time scales are of the same order. Employing a master equation approach, 
we obtain the stationary reduced density matrix retaining the coherences between vibrational states. 
Depending on the system parameters, two regimes are identified, characterized by either (i) an 
effective thermal state with a temperature lower than that of the environment or (ii) strong coherent 
effects. A marked cooling of the vibrational degree of freedom is observed with a suppression of the 
vibron Fano factor down to sub-Poissonian values and a reduction of the position and momentum 
quadratures. 

PACS numbers: 85.85. +j, 73.63.-b 



I. INTRODUCTION 

Nano-electromechanical system^ (NEMS) represent 
an intriguing class of devices composed of a nano- 
mechanical resonator coupled to an electronic nanode- 
vice. Several examples of NEMS have been realized, 
ranging from single oscillating molecules,^ to suspended 
carbon nanotubeJ^^'^ and suspended nano-cantilevers or 
nano-beams.'^ 

In all these devices, the coupling between mechanical 
and electronic degrees of freedom gives rise to peculiar 
transport phenomena like Franck-Condon blockade,^'^^!^ 
negative differential conductance^^^^M^ and remarkable 
noise characteristics.^^ Along with their outstanding 
electronic properties, also mechanical ones are of extreme 
interest. Indeed, owing to the extreme sensitivity of the 
vibrating part to the spatial motion, NEMS have been 
proposed as novel detectors in scanning microscopeJ^^ 
or as ultra-sensitive nano-scales, being able to detect 
the mass of even few molecules adhering to them.l^^I^ 
In order to successfully employ a NEMS as a precision 
position detector it is important to reduce its thermal 
fluctuations, eventually attaining the ultimate goal of 
cooling it down to its quantum ground state.^^ Also, 
ultra-sensitive NEMS position detectors based on pecu- 
liar quantum states such as position-squeezed state^^ 
have been proposed and experimentally realized. 
A great variety of NEMS setups have been investigated 
theoretically. Typically, the electronic part is composed 
of a semiconducting normaP^^^ or supercon- 
ducting^^ single quantu m dot .^^ Double-dot setups 
have been studied as well.^^'^ Different models for 
the coupling between electrons and vibrons have been 
considered, ranging from the simple Anderson-Holstein 
(AH) model^ to microscopic models tai- 

lored for specific systems, such as suspended carbon 
nanotub es. '^^ "^^ ' Also the influence of external dissipative 
baths^^ 46|47l radiation fields^^ has been analyzed. 
In the most complex configurations, interesting physical 
effects have been predicted. For instance, for a nano- 



mechanical resonator coupled to a microwave cavity,'^ 
to field driven quantu m do^ ^ or in configurations with 
double quantum dots,'^^'^^' phonon cooling has been 
found. With radio frequency quantum dots^^ and a elec- 
tromagnetic cavit}'^^ squeezing of the vibron position and 
momentum quadratures has been theoretically predicted. 

Even the simple AH model for a single electronic level 
coupled to an undamped vibrational mode exhibits a rich 
physics, part of which is still unexplored. Roughly speak- 
ing, two regimes have been considered so far, according 
to the ratio 
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between the vibron frequency, cjq, and the average dot- 
leads electron tunneling rate, Fq. 

For fast vibrations, 7 <C 1, every electron tunneling event 
occurs over many oscillator periods. Then the electrons 
are not sensitive to the position of the oscillator, but 
only to its energy.^^ Consequently, the oscillator den- 
sity matrix becomes close to diagonal in the basis of the 
energy eigenstates.^^ Many groups focused their atten- 
tion on this regime, employing rate equations to study 
transport phenomena as th e Franck -Condon blockade, 
super- Poissonian shot noise I 9 | 1Q | 12 | 24 | even peculiar 
effects such as sub-Poissonian out of equilibrium vibron 

distributions.28 29 33 56 57 

In the opposite regime of slow vibrations 7 ^ 1, 
electrons are extremely sensitive to the position of 
the oscillator, which can be treated in a semi- classical 
approximation.^^'^^'^^ll] indeed, Mozyrsky et al. have 
shown^^ that it is the onset of a semi-classical Langevin 
dynamics. In this regime the electronic properties of the 
system have been especially investigated, in particular 
the current and shot noise^^ in both limits of wealP^ 
and strong electron- vibron coupling.'^ In the latt er case, 
bi-stability and switching have been addressed. SSl^IIlol 
Recently, the classical phase space of the vibron has also 
been studied.^ 
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Less attention has been devoted so far to the coherent 
regime, where the off-diagonal elements of the system 
density matrix in the energy representation play a rel- 
evant role. In this regime, which star ts arou nd 7^1 
(for a more precise discussion, see Sec. Ill A), the com- 
petition between the vibron and the electron time scales 
gives rise to a tough theoretical problem. Most of the 
results obtained so far, concerning bi-stability and phase 
space analysis, have been obtained stretching somehow 
the validity range of semi-classical approaches j ^^ * ^^ * ^ -^ in 
the limit of very low temperatures, r <C 1, where 



kBT 



(2) 



here is the Boltzmann constant and T the environ- 
ment temperature. The interplay of electron and vibron 
time scales is expected to strongly influence the dynam- 
ics of the vibron. For instance, for a NEMS based on a 
metallic dot in the weak coupling regime, the damping 
effect of tunneling electrons on the vibron dynamics is 
maximal when ujq ^ Fq.^ Similar mechanisms could 
play a role also in the simpler model of a single level 
quantum dot. 

Motivated by these considerations, in this article we 
investigate the vibronic properties in the coherent regime 
7 > 1. Here, because of the off-diagonal structure of 
the system density matrix, a simple rate equation is no 
lon ger^ ju stified.^^ We derive a generalized master equa- 
tion^^^^ in the sequential tunneling regime, in which 
all off-diagonal elements of the reduced density matrix 
in the energy eigenbasis are retained. In the limit of 
high temperatures considered here, r > 1, a fairly large 
number of basis states have to be included. This fact 
leads to a serious numerical challenge. Our calculation 
extends up to r < Tmax, where Tmax ^10. We remark 
that, as a difference with previous studies our ap- 
proach is not restricted to small electron- vibron coupling. 

Here is a summary of our findings. With the ex- 
ception of a region 7 ^ r, in the stationary regime the 
vibron state can be approximately described in terms of 
an effective thermal distribution. In the coherent regime, 
the effective temperature is lower than the environmen- 
tal temperature. Here, the role of coherences is crucial 
despite subtle. Non- vanishing off-diagonal elements of 
the vibron density matrix in the eigenbasis, despite 
being very small compared with diagonal elements, 
originate the peculiar effective thermal re-distributions 
of the vibron occupation probabilities. 
When ^ the system exhibits deviations from the 

above effective thermal state. The off-diagonal elements 
are larger, leading to a marked suppression of the 
vibron fluctuations even below the Poissonian value. 
The main results of our paper concern the stationary 
vibron properties in the coherent regime and can be 
summarized as follows: 

(i) a cooling of the vibrational mode with respect to the 



temperature of the electronic environment; 

{a) a strong suppression of the vibron Fano factor 

eventually reaching sub-Poissonian values; 

(Hi) a reduction of the variances of the vibron position 

and momentum quadratures. 

We remark that the reported cooling phenomenon is a 
direct consequence of the NEMS entering the coherent 
regime. It is not "induced" by any external drive 
or dynamics, like connecting the system to several 
reservoirs at different temperatures.^Sl 
All the above effects are more pronounced when the 
electron-vibron coupling strength is increased and none 
of them comes out treating the system with a simple 
rate equation involving the diagonal matrix elements 
only. 

The paper is structured as follows. In Sec. [H] we describe 
the Anderson-Holstein model and the derivation of the 
generalized master equation in the stationary regime. In 
Sec, mil we illustrate the coherence effects on the vibron 
behavior and on the electronic degree of freedom. The 
large discrepancy with respect to the results obtained by 
means of a simple rate equation is highlighted. Conclu- 
sions are drawn in Sec IIVI 



II. MODEL AND METHODS 



Anderson-Holstein Model 



In the AH mo del, ^I^'^^'^ the hamiltonian 



H = K 



dot 



(3) 



describes a ultra-small quantum dot (i^dot) coupled to 
an harmonic oscillator (Hpsc) yia- the coupling term i^int • 
The quantum dot is modeled^ as a spin degenerate single 
level with the average level spacing of the order of the 
charging energy Eq (from now on, h = 1) 



Hdot = eh^Ech{h- 1). 



(4) 



Here, n = X]cr=±i is the occupation number of the 
level, with her = d^^d^ the occupation of spin a/2 with 
cr = ±1 and dcr, d}^ are the fermionic dot operators. We 
assume that Eq is the largest energy scale of the problem 
and we consider only single excess occupancy on the dot, 
n = 0, 1. The energy e = ^ + 2Ec{l/2 — rig) includes the 
energy of the lowest unoccupied single-particle level ^ and 
a term connected to ng = CgFg/e, the charge induced by 
the gate voltage with gate capacitance Cg (— e is the 
electron charge). 

The vibron is described as an harmonic oscillator with 
mass m and frequency ujq. In terms of the boson opera- 
tors 6, it is modeled as 



^01 



(5) 
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The dot and the oscihator are coupled via a term bi-hnear 
in the oscillator position x and in the effective charge 
number on the dot, h — n^^^^ 



Hint = V2Xujo^{n - rig) , 

^0 

where A is the adimensional coupling parameter and 

1 



(6) 



(7) 



is the characteristic length of the harmonic oscillator. 
The dot is coupled also to the external left (L) and right 
(R) leads of non interacting electrons 

Pleads =^ ^ ^a^k,a^c^,k,a , (8) 

a=L,Rk,a=±l 

where Ca,k,a and cj^ ^ ^ are the fermionic operators. The 
leads are assumed in equilibrium with respect to their 
electrochemical potential /iL,R = Mo i where V 

is a symmetrically applied bias voltage, and jno is the 
reference chemical potential. The bias V forces electrons 
to flow from the left to the right lead through the dot via 
the tunneling hamiltonian 



= ^oE E 



^a,k,a^o 



h.c. , 



(9) 



a=L,R/c,cr=±l 



where to is the tunneling amplitude through both the left 
and right barriers. 

Equation (|3| can be diagonalized by the Lang-Firsov po- 
laron transformation,'^ with generator 

U = exjp [fj {b^ — b)] and f} = \{fi — ng). (10) 

This procedure imposes no restriction on the possible val- 
ues of A. The transformed operators in the polaron frame 
O = UOW SiTeb = b-fi and 4 = exp [X{b - b^)] , 
while h is invariant. The diagonal Hamiltonian, ex- 
pressed in terms of the original operators reads 



H 



(11) 



with renormalized level position e = ^ + {Eq — A^cJo)(l — 
2ng). In the following we choose /io = ^ setting the 
resonance between the n = 0, 1 states at rig = 1/2. The 



eigenstates of Eq. (11) will be denoted as |n,/), where n 
is the dot occupation number and / represents the vibron 
number. The transformed tunneling Hamiltonian is 



Hi — to 



E 



E 



A(6-fct) t ^ 



h.c. 



(12) 



a=L,Rk,a=±l 



B. Master Equation 

The dynamics of the dot and the oscillator is described 
by the reduced density matrix deflned as the trace 



over the leads of the total density matrix ptotC^), in the 
polaron frame 



p{t) = Trieads{Ptot(^)} . 



(13) 



We consider the sequential, weak tunneling regime, treat- 
ing Ht in Eq. (12) to the lowest order. This approxi- 



mation is valid for not too low temperatures T, Fq = 
27r|topi^ < kBT^ where u is the leads density of states. 
We further perform the Born approximation,^^ assum- 
ing that the system and the leads are independent be- 
fore Ht is switched on. This amounts to take the fac- 
torized form at t = 0, Ptot(O) = p(0) ^ Pi{^)^ where 
pi(0) = Pl(0) (8) Pr(0) and Pl/r(0) are the initial equi- 
librium density matrices of the left and the right lead 
respectively. 

In the interaction picture with respect to Ht^ any op- 
erator, A, is transformed as 

The reduced density matrix evolves according to 

Piit)=-J2 /dt'|[Qi..(t),QiV(*')Pi(*')]^+ 

-[QiAthPi{t')QUt')]K-{t' -t) (14) 
+ [Qljt),Qi,„{t')pi{t')]K-{t-t') 
- [QlAt),p,{t')Q,At')]K+{t'-t)} . 

Here, 

and K^{t) = K^{t) + K^{t) are the leads correlation 
functions 

Kit) = ItoPE^rieads {clkA^)Ca,kA^)pi} , 
k 

K-{t) = |to|'^TVieads{ca,fe,.W4^_^(0)pi} . 



(15) 



In Eq. ( 14 ) we performed the large reservoirs approxima- 
tion^ 



Pi,tot(^0 = Pi{t') • Pi 



(16) 



assuming that tunneling events have a negligible effect on 
the leads, which remain in the thermal equilibrium state, 
denoted as p\. In the weak tunneling regime (Fq < k^T) 
one can replace pi(t') ^ pi(t) using the standard Markov 
approximation, and extend the integration limit to oo. 

Eq. (14) can be projected on the eingenstates of the 
hamiltonian (pTl), obtaining an inflnite set of coupled 



equations for the density matrix elements (n, l\p{t)\n' ^ V) 
(where n, n' G {0, 1} and /, T > 0). It can be easily shown 
that diagonal and off-diagonal elements in the electron 
number decouple and in the stationary regime the latter 
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tend to zero. In fact, the coupling of the leads and the 
dot charge leads to a rapid d ecay of superpositions of dot 
states with different charges !^^^ Since we are interested 
in the stationary properties, we focus on the density ma- 
trix elements diagonal in the electron level occupation 
number, p^qf{t) = {n^q\p{t)\n^q') . They obey the follow- 
ing generalized master equation (GME) 

K,At) = E 7^^,Ve"^"'''"'''"'^''^* Ptw^t) ' (17) 

pp'n' 

where TVl'^,^^, are the Redfield tensor element 



-\-XiqXipC {UplYSj 



^Iq'pp' —'^XpqXpfqf [C+ (cjgp) + C+ (cjg/p/ ) *] 



i^qqippi 



Here 



^-2Y,{Xq,lXp,lC^{uJlp.y^p 



-VXqiXpiC^ {{xllp)bpi q^ 

■ uJo{q — q') and 



(18) 



(19) 



(20) 



is a Fermi function, (E) = 1 — f^{E) and 



^a{^qq') = log ( ^ ) + Rc 



1 iP ^_ 

2 + + 



Ma) 



the depolarization shift. The stationary system proper- 
ties are obtained from the solution of the steady-state 
GME 



^ppf 



0, 



(24) 



pp'n' 



written here in the Schrodinger representation, for the 
stationary values of the reduced density matrix elements 



(25) 



Note that the GME p4| ) takes into account both all off- 
diagonal (coherences) and diagonal terms in the vibron 
number. 

The steady state current calculated e.g. on the right 
tunnel barrier, reads 

I = eTo ^ {-2 [C+M + C+M*] X,pX,,pp°^, 

pqq' 

+ [CnM + C^M*] Xp,Xp,,pl^,} .(26) 

Note that in the steady-state the current is independent 
of the barrier index. 



are the generalized Franck- Condon factors^^, with q^ = 
mm {q^q'}^ q^ = mdix{q^q'} and Lq{x) the generalized 
Laguerre polynomials.^^ The factors 2 in Eqs. (19) are 



due to the spin degeneracy. In Eqs. (18),(19) the gener- 
alized tunneling rates {uqq' ) = C^l^q^^ {uqq' ) 
have also been introduced 



POO 

Jo 



(21) 



Exploiting the identity 

) 

d6'exp {ine) = n5{n) + iP.V.{l/n) , 



/ 

Jo 



and the explicit form of the leads correlation functioiP^ 



with Uc the cut-off energy and f3 = 1/(/cbT), one has 

a. 

Here 



1 



1 ^ ^(3{e^E-Hc.) 



(23) 



C. Rotating Wave Approximation 

In the regime of fast vibrational motion, 7^1, the 
contribution of fast oscillatory terms to the solution of 
Eq. (17) is negligible. It is then possible to perform the 



rotating wave approximation (RWA), neglecting the os- 
cillating contributions and including only the dominant 
secular terms. They are those which connect density 
matrix elements with p—p' = q — q'- This implies that the 
diagonal elements pqq are decoupled from the off-diagonal 
ones Ppp, with p^ ^ p. In addition, the non-diagonal ele- 
ments vanish in the stationary regime. Hence station- 
ary properties are fully described by the diagonal occu- 
pation probabilities Pnq = pqq and Eq. (24) reduces to a 
standard rate equation 



ZnPnq ^ ^gp 



EE 

n'y^np=0 n'y^np=0 



2„,P„,pr""=0. (27) 



The coefficients Zn stem from the spin degeneracy: zq = 
2, zi = 1. The tunneling rates T^^' = Tl^^ + T^^^^^ for 
the transition |n,p) |^^<7), are 

r^^gp^ = 2XpqRe {C^ {ujpq)} = ToXpqf^ {ujpq); ^^^^ 

^aVq^ = '^XpqRe {C~ {ujpq) } = V^X^qf' {ojpq). 
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FIG. 1: (Color online) Vibron occupation probabilities Pq as a 
function of the vibron number calculated numerically solv- 
ing the RWA rate equation in Eq. ( [27| (red empty squares) 
and the GME in Eq. (24). Here 7 = 3 (cyan empty tri- 
angles), 7 = 1 (purple empty circles), 7 = 10~^ (blue cir- 
cles), 7 = 10~^ (green squares), (a) A = 0.5 and eV = cjq- 
(b) Main panel: A = 0.5 and eV = 20ujo. Inset: A = 2 
and eV = 20(jJo. In all panels other parameters are r = 6, 
rig — 1/2 and Uc — IO^cjq- 



Note that Eq. (27) does not depend on 7. 



Within the RWA, the steady-state current Eq. (26) is 



^RWA 



0,1 



(29) 



Results presented in the following Section are obtained 
by numerical solution of the GME Eq. ([24| truncating 
the size of the vibron Hilbert space until convergence 
is reached (see next Section for details). In the limit 
7 <C 1, this approach accurately reprod uces the w ell- 
known solution of the RWA rate equation.l^ ^^ ' ^^ ' ^" ' 

This is illustrated in Fig. [T] where the diagonal elements 
of the reduced density matrix, Pq = p^^, obtained by 
solving the GME for different values of 7 are shown. The 
solutions of the GME converge to those in the RWA for 
7^0 (red squares), both for small and large voltages 
(panels (a) - (b)) and even in the strong electron- vibron 
coupling regime (panel c). This convergence has been 
systematically observed in the whole range of parameters 
and constitutes a validation of our numerical procedure. 
On the other hand, with increasing 7, deviations from 
the RWA are obtained. They are originated from the 
coherent off-diagonal elements of pqq>. These deviations 
represent the central part of our work and relevant conse- 
quences will be discussed in details in the following Sec- 
tion. 




FIG. 2: (Color online) (a) Stability diagram in the (1/, rig) 
plane: gray shaded regions denotes Coulomb Blockade regime. 
Blue (solid) lines mark the onset of transitions between 
ground states, red (dashed) lines signal transitions involving 
excited states of the vibron. (b) Different regimes as a func- 
tion of 7 = Fo/cJo and r = /cbT/cjq- The blue shaded region 
denotes the "diabatic" regime of fast vibrations, the red one 
the semi-classical regime with cjq ^ Fq. The region above the 
line 7 = T marks the region accessible by our methods and 
the yellow shaded area indicates the accessible regime where 
coherences among vibron states become important. 



III. RESULTS 

In the present Section, we investigate the stationary 
equilibrium properties of the system as a function of the 
parameters 7, r defined in Eqs.([T]), ([2|. Here, 7 distin- 
guishes between fast, slow or coherent regimes, and r is 
the reduced temperature. We focus mainly on the vi- 
bronic properties. Relevant quantities of our interest are 
the vibron Fano factor and the position and momentum 
quadratures. Their behavior can be understood by first 
analyzing the structure of the oscillator density matrix 
in the vibron eigenbasis. In the final part of this Section, 
the electronic current and the average electronic occupa- 
tion of the dot will be briefly addressed. 



A. Parameters regimes 

Figure [2|a) represents the stability diagram of the 
system as a function of the bias voltage V and of rig. In 
the shaded (gray) regions the system is in the Coulomb 
blockade regime with the dot empty (full) if ng < 1/2 
(ng > 1/2). The blue (solid) lines mark transitions 
between the states |0, q) and |1, q) {n and q are the occu- 
pation number of the electronic or vibronic states \n^q)). 
Red (dashed) lines indicate the activation threshold 
for transport channels with transitions \n,q) \n'^q') 
where in general q.,q' 7^ 0. The boundaries of different 
transport regions are thermally broadened. 
In the following analysis we will focus at ng = 1/2 where 
the n = and n = 1 charge states are on resonance 
(dashed-dotted line in Fig. [2|a)). Qualitatively similar 
results are observed also off-resonance at rig ^ 1/2. 
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In Fig. [2|b) we report a sketch of the system behavior 
in the (7, r) plane. 

We distinguish a ^^diabatic^^ regime of fast vibrations 
for 7 <C 1 (blu e shaded region), where the RWA is 
^^l^^mmEmm and an ''adiabatid' regime of slow 
vibrations 7 ^ 1 w here semiclassical methods have 
been applied^^H^^'^^'^ with approaches confined to low 
temperatures r <C 1. 

The regime in between these two regions is, on the other 
hand, unexplored. Our GME method fills in precisely 
this gap. Indeed, we will show that, increasing 7 from 
the diabatic regime, marked deviations from the solution 
of the RWA rate equations appear. Here, coherences 
represented by the off-diagonal elements of the reduced 
density matrix are relevant. The constraints on our 
GME method are: 

{i) the sequential tunneling approximation which implies 

< 7 < r (the area above the red line in Fig. |2|; 

(ii) a restriction on the temperature r < Tmax ^ 10, due 

to the rapid increase with r of the number of vibron 

states needed for the convergence of our numerical 

calculations. 

Because of these constraints, the transition towards the 
semiclassical limit cannot be explored. 

We solve numerically the GME in the steady state, in- 
creasing the size of the vibron Hilbert space including 
up to 150 oscillator states per charge state in the ba- 
sis, until convergence is reached. The limitation on the 
number of oscillator states employed in the calculation is 
essentially due to computer memory requirements and by 
the decreasing rate of convergence. At high temperature 
r ~ 10, already more than 100 oscillator basis states are 
required for convergence. 



B. Overview of the results 

Before entering our analysis, we summarise in Fig. |3] 
the stationary vibronic characteristics we obtain in the 
various parameter regimes. We address the temperature 
regime 1 ^ r < Tmax and < 7 < r. We will show that in 
a large parameters range (shaded area) the vibron density 
matrix 



Pqq' 



n=0,l 



(30) 



is well approximated by an effective thermal distribution 
with temperature 



-(th) / X 



(1 _ e-i/-e„ j 



(31) 



For 7 <C 1 the off-diagonal elements of the reduced den- 
sity matrix are negligible and a fit of pqq^ to Eq. (31) 
leads to Teff > r. This "heating" phenomenon is due 
to the finite voltage, and for V ^ the system attains 
a thermal equilibrium distribution in the polaron frame. 




FIG. 3: (Color online) Schematic overview of the results in 
the (7, r) plane, see text for details on the notation. 



'Teff ^ r, as reported in Ref. 

With increasing 7 the off-diagonal elements of pqq' in- 
crease. The onset of this coherent regime is marked by 
the condition 7 > 7(A), with the latter a A-dependent 
threshold value. For typical values A ^ 1 we find 
0.1 < 7(A) < 1. When 7(A) ^ 7 < t the coherences 
are much smaller than the diagonal elements of the re- 
duced density matrix. In this weakly coherent regime the 
vibron density matrix can still be fitted with a thermal 
distribution, but at a lower effective temperature, even- 
tually reaching the cooling regime where Teff < r. The 
cooling is always accompanied by a reduction both of the 
fiuctuations of the vibronic population and of the vari- 
ance of position and momentum quadratures. 
A completely different system behavior takes place when 
7 ^ r, (yellow-green area delimited by the dashed 
and continuous lines in Fig. [3|. In this strongly coher- 
ent regime off-diagonal terms of the density matrix are 
paramount and the diagonal part of the density matrix 
deviates from the simple thermal distribution. Here, de- 
spite of the high temperature (r > 1), a non-classical 
system behavior comes out. Suppression of the vibron 
populations fiuctuations and of the variances of the po- 
sition and momentum quadratures persists also in this 
regime. 



C. Density matrix and cooling 

We start our analysis investigating the structure of the 
oscillator density matrix in the vibron eigenbasis in the 
regimes indicated in Fig. [3| 

We first concentrate on the density matrix pqqf when it 
is well described by the effective thermal distribution, 
Eq. ( [31] ). In Figs. [4]^ a) - (b) we report the density plots 
of \pq^ in the diabatic and in the coherent regime re- 
spectively. For 7 <C 1 - panel (a) - the density matrix is 
strongly peaked around the diagonal, q ^ q' . The cor- 
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FIG. 4: (Color online) Density plot of \pqq' \ as a function of 
q^q'. Panel (a): 7 = 0.01 and r = 9; Panel (b): 7 = 2 and 
T = 9. Panel (c): Occupation probabilities Pq = pqq for r = 9 
and 7 = 0.01 (squares) or 7 = 2 (dots). Lines are fit to a 
thermal distribution with effective temperature Teff. In panel 
(d) 7 = 1.2 and r = 3 (main), and r = 9 (inset). In all panels, 
rig — 1/2, eV = cjo, A = 2 and uoc — 10^ cjq- 



responding occupation probability distribution Pq = pqq 
extends over more than fifteen states Fig.[4]^c) (squares). 
We then perform a numerical fit of pqqf on a thermal 
distribution in Eq. (31). The fit leads to an effective 



temperature Te// with an error Ar, signaling the depar- 
ture from the approximatively diagonal thermal density 
matrix. In considered diabatic regime, 7 <C 1, Teff ~ r 
with a very small relative error 5r = Ar/reff < 10~^. 
Increasing 7, the coherent regime is entered, with a 
vibron occupation probability considerably altered, see 
Fig.[4]^b). We observe two main modifications. First, off- 
diagonal elements are larger than in the diabatic regime, 
even though they remain rather small in comparison with 
the diagonal ones. Second, the probability distribution 
along the diagonal gets narrower. Remarkably, the occu- 
pation probabilities Pq are still approximated by a quasi- 
thermal state but with an effective temperature Teff lower 
than the environmental one, see Fig.jlj^c) circles. The rel- 
ative error is still rather small, ~ 6 • 10~^. We remark 
that this result is a consequence of the non- vanishing vi- 
bronic coherences, in fact in the RWA the density matrix 
would be strictly diagonal, as a difference with Fig. [Ij^b). 

The crossover from heating to cooling with increasing 
7 and the accuracy of the thermal approximation mea- 
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FIG. 5: (Color online) Solid lines: effective temperature Teff 
as a function of 7, extracted from a numerical fit of pqq^ on 
Eq. (|3l]), (see text). Panel (a): r = 9 and eV = 20cjo. In 
the inset Teff as a function of A for 7 = 6. Panel (b): r = 3 
and eV = Gcjq- In all panels, dashed lines around the solid 
line delimit the absolute error on the effective temperature, 
Ar. Red (blue) gradient areas signal heating (cooling). Other 
parameters are rig = 1/2 and cjc = 10^ uoq. 



sured by the error of the fitting procedure are analyzed 
in Fig. [sja) - (b). The absolute error is represented by 
the shaded area around the continuous line limited by 
dashed lines. Two different regimes are clearly identified. 
In the diabatic regime, the effective temperature is larger 
than that of the environment (red shaded region) and 
the system exhibits heating. The values of Teff > r are 
due to the considered high voltage bias and indeed for 
V ^ one obtains Teff r. 

As 7 is increased, cooling occurs (blue shaded region) 
and Teff drops markedly below the electronic tempera- 
ture r. Even though Teff increases for increasing V (not 
shown) this cooling effect survives up to eV ^ 20cc;o. 
It can be seen that the error grows with increasing 
7, signaling the rise of the off-diagonal terms of the 
density matrix. We can then safely identify the cooling 
phenomenon provided 7 does not approach r. In the 
inset of Figjsja) the typical behavior of Teff as a function 
of A is shown, with an error on Teff which increases 
increasing A. 

We remark that the cooling phenomenon is entirely due 
to the coherent dynamics of the vibron-electronic system 
and is not induced by any ad-hoc mechanism acting on 
the system. 

The description in terms of an effective thermal distri- 
bution ceases to be valid when 7 ^ r, as shown by the 
increasing errors in Fig. [5] To investigate this regime, 
we consider the case r = 3, 7 = 1.2 in Fig.[4];d). Here, 
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although an effective temperature Teff ~ 0.8 < r can be 
formany extracted with the fitting procedure, the rela- 
tive error becomes large 5t ~ 0.11. The main source of 
error are the rather large off-diagonal matrix elements. 
This is a general trend which we always observed when 7 
approaches r. As an illustration, we report in the inset 
of Fig. [4]^c) the density matrix for the same parameters 
of the main panel but at higher temperature (r ^ 7). In 
this case off-diagonal elements are strongly suppressed 
and an effective thermal description is then appropriate. 
We conclude this paragraph commenting on the depen- 
dence of the above results on the electron-vibron coupling 
strength. A, which is not limited in our approach. In Fig- 
urel6](a) - (d) we report the density plot of \pqq> \ at fixed r 
and for increasing values of A. Increasing the dot-vibron 
coupling induces an hybridization of the electronic and 
mechanical degrees of freedom which manifests itself also 
in the off-diagonal elements of the vibron density matrix 
between Fock states. Indeed, with increasing A a "de- 
localization" phenomenon is induced, the density matrix 
spreading away from the diagonal. 



tioiPl 







\Pqq' 



0.1 



Hi (b) 



(c) 








FIG. 6: (Color online) Density plot of \pqq' \ as a function of 
q for fixed 7, r and different values of the electron-vibron 
coupling strength. Panel (a): A = 1; Panel (b): A = 1.5; 
Panel (c): A = 2; Panel (d): A = 2.5. In all panels, rig = 1/2, 
eV = cjo, 7 = 3, r = 6 and ujc — 10^ cjq- 



D. Wigner function and quadratures 

Further insight on the system behavior is obtained 
from the Wigner quasi-probability distribution func- 



1 

W{x,p) = - / d^ e^'Py{x - y\U^pU\x + y) , (32) 
^ J -00 

the quantum analogue of the Liouville density in classi- 
cal phase space with U defined in Eq. ([lO| . The Wigner 
function allows the detection of non-classical features sig- 
naled by regions in the {x^p) plane where W{x^p) < 0.^ 



0.02 




W{x,p) 



0.1 




P{x) 



FIG. 7: (Color online) Panels (a) - (c) density plot of the 
Wigner function W{x^p) as a function of In panels (a) 
and (b) r = 9, and in (a) 7 = 0.01, in (b) 7 = 2. Note 
the different color scale in the two panels, (c) Density plot of 
W{x^p) for r = 3 and 7 = 1.15. (d) Probability density of the 
oscillator position P(x). In all panels x and £^are expressed 



in unity of Iq and Iq respectively - see Eq. (|7|. In addition 



1/2, = cjo, A = 2 and uJc 10^ 



Figures [t]^ a) - (b) show W{x^p) for the parameters of 
Figs. |4]^a) - (b), in the effective thermal regime. The 
Wigner functions are positive and can be regarded as the 
probability density of the oscillator states in the (x,p) 
phase space. Consistently, their shape closely resembles 
that of a thermal state, but with a decreased width 
Teff. The shrinkage of W{x^p) observed when entering 
the coherent regime can be traced back to the role of co- 
herences which mediate the redistribution of the Pq and 
the ensuing reduction of Teff. 

Figure [tJ^c) shows W{x^p) for 7 = 1.2 and r = 3 in the 
strongly coherent regime, when the system can no longer 
be described by a thermal distribution (as in Fig. [4jd)). 
In this case, the width of Wix^p) is still reduced. As 
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a difference with previous case, here we observe that 
W{x^p) becomes negative in an approximately circular 
area around the main peak (white dashed circle). This 
fact signals a non- classical hehdivioi of the system entirely 
due to the non-diagonal structure of the reduced density 
matrix. Negative values of W{x^p) are a trademark of 
the regime 7/r ~ 1. This quantum behavior, naively un- 
expected at the considered high temperatures r > 1, is 
a relevant result whose implications on the vibronic fluc- 
tuations will be discussed in the next subsection. 
Also noteworthy is the elongated shape of W{x^p) along 
the X direction, Fig.[7|^c). This is reflected in the position 
probability distribution 



P{x) 



dp W{x^p) 



(33) 



shown in Figure ^d). While in the effective-thermal 
regime P{x) exhibits a single peak, in the regime of 
strong coherences it shows a peculiar shoulder structure. 
In this latter regime, for relatively slow tunneling events 
7 ~ 1, the oscillator adjusts itself to the two equilibrium 
positions corresponding to zero or one extra electron on 
the dot, xq = Xug and xi = A(n^ — 1) (in units of io^ see 
Eq. When the width of the two probability distri- 

butions for the uncharged and charged states (ex Teff) is 
of the order of the separation between the rest positions 
(X A, the shoulder is visible. The effect is therefore more 
pronounced the larger is the electron-vibron coupling A 
and eventually develops into a double-peaked shape. 
The width of W{x,p) in the x and the p directions is 
strictly related to the variance of the vibron position and 
momentum 



dx 



dp x^W{x,p) , 



(34) 



and corresponding expressions for p. In Fig[8](a) and (b) 
are reported the variance of the position, X = x/£o = 
{b^ + b)/V2, and momentum V = p • io = - b)/V2 
satisfying always Var(A')Var(7^) > 1/2 - with £q defined 
in Eq. As expected, in the coherent regime both vari- 
ances are strongly suppressed with respect to the large 
values attained in the diabatic regime, both for high and 
low (not shown) voltages. This fact denotes a tendency 
towards squeezing^ 



E. Vibronic Fluctuations 

The out-of-equilibrium vibron behavior is conveniently 
discussed introducing the vibron Fano factoi'^^'^ 



Var(7V) 



(35) 



where N = b'^b. Here, is the ratio between the 
variance of the vibron occupation number Var(A/') = 
(Ar2) - (Ar)2 and its average {N), where {O) = Tt{Op} = 
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A = 2 
■A = 0.5 



0.01 
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FIG. 8: (Color online) Variance of the position (a) and of 
the momentum V (b) as a function of 7. In all panels r — 9, 
eV = 20a;o, '/T'g = 1/2 and uoc — 10^ ujq. 



Tt{UOW p}. This quantity (not be confused with the 
electronic Fano factor) brings information about the 
statistics of the vibronic mode and also on electronic 
properties, like charge fluctuations. This is due to the 
involved polaron transformation and it can be explicitly 
seen introducing the hybrid average {0)p = Tr{Op}. In 
terms of {■)p one has 

(V)=(V),--V2(7)A') + (f) 
Var(V) = V{N) + V{f]^) + 2V(7)A') + 2C(7)^ V) (36) 
- V2 [C{fjX, N) + C(V, fjX) + 2C(f , fjX)] , 



where V(0) 



{0)l and C{A,B) 



{AB)p — {A)p{B)p, with adimensional oscillator po- 
sition X. The quantities (N) and Var(V) depend 
explicitly both on the charge (being fj = X{h — Ug)) 
and on its fluctuations. In Eq. ([36| all terms of the 
form {■ X) and C(-, •A') depend on the coherences of 
the density matrix. In the RWA these terms simplify 
but keep the 7) dependence. We remark that, beside 
the explicit dependence of Var(V) and {N) on (n), 
fluctuations of the electronic charge intrinsically affect 
Fv, since the contribution of electronic and vibronic 
degrees of freedom cannot be factorized in p^^,. 
Fv belongs to the class of "bosonic" Fano factors, 
originally introduced to characterize boson distributions 
and often used to study photon populations in the 
context of quantum optics. Typically, for bosons the 
Fano factor is super- Poissonian {F^ > 1). This is, 
for example, the situation for light from a classical 
radiation field.'^ Deviations from the super-Poissonian 
regime signal non-classical correlations, which have been 
observed for instance in experiments with micromaser.l^ 
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From an experimental point of view, the situation 
for vibron is more complicated than in the case of 
photons, although many proposals have been made to 
observe quantities like the average vibron numberP^H^ 
Recent experiments aiming at characterizing the vibron 
distributions have appeared!^ A few recent theoretical 
works investigated the vibron distributions predicting 
sub-Poissonian vibron Fano factor in the diabatic regime 
7 <C 1, for low temperatures r <C 1 and in selected 
parameter ranges. 
The behavior of the vibron Fano factor as a function 



quantum behaviors show up in the steady state regime. 
The vibron Fano factor also depends on the voltage and 
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FIG. 10: (Color online) Main: plot of as a function of V 
(units cjo/e) for different values of A. Inset: dependence of 
Fv on the electron- vibron strength A, for eV = IScJo- Other 
parameters are 7 = 6, rig = 1/2, r = 9 and ujc — 10^ cjq- 



on electron- vibron coupling. increases with increasing 
bias voltage, assuming super-Poissonian values even in 
the coherent regime. Fig. [To] (main panel). This behavior 
can be attributed to the increase of scattering between 
oscillator states induced by the current flow across the 
system. In the inset, the Fano factor is plotted as a 
function of A: decreases with increasing coupling 
strengths. This trend has been found both in the low 
{eV < k^T) and in the high (eV > k^T) voltage regimes 
- see also Figs, ^a) - (b). 



FIG. 9: (Color online) Vibron Fano factor F^ as a function of 
7 for different values of the electron- vibron coupling strength 
A and temperature r. In panel (a) r = 9 and eV = cjq, in 
the inset eV — 20 cjq- In panel (b) r = 3 and eV — ujq. 
Solid lines are the results obtained with the GME, dashed 
lines are obtained within the RWA. In all panels, rig = 1/2 
and uJc — 10^ ujq. 



of 7 for weak and strong A is reported in Fig. [9j In 
the regimes where the vibron can be approximated 
by an effective thermal distribution, F^ displays a 
super-Poissonian behavior. This occurs in the diabatic 
regime, both at low and high voltages (eV smaller or 
larger than k^T)^ and in the weakly coherent regime. In 
this case however the vibron Fano factor is considerably 
reduced. These behaviors qualitatively do not depend on 
the electron-vibron coupling A. In the diabatic regime 
the RWA applies with no dependence on 7. 
On the other hand, entering the strongly coherent 
regime (7 ^ 1 and r = 3 in Fig. [ofb)), non-classical 
sub-Poissonian vibron Fano factors are obtained. This is 
in correspondence with the negative values of the Wigner 
function observed in this parameter regime. Both fea- 
tures result from the large vibronic coherences. Similar 
conclusions were drawn from the frequency-dependent 
shot noise of a NEMS.^^ Here, however, these peculiar 



F. Charge degree of freedom 

We conclude the survey of our results briefly com- 
menting on the electronic properties. Figure [TTJa) shows 
the I(y) characteristics for different values of 7. In the 
regime r > 1, the vibron sideband can not be resolved 
and the current increases monotonically. Figure pTj^b) 
shows the ratio between the current / obtained from the 
GME and the current in the RWA, /rwa, as a function of 
7 and fixed voltage V . Increasing 7, the current increases 
with respect to the RWA limit. This qualitative trend oc- 
curs in all the parameter regimes explored. The increase 
of the current is more prominent for larger values of A. 
This is a consequence of the delocalization of the vibron 
density matrix occurring with increasing electron-vibron 
coupling. 

The average occupation of the electronic dot level (n) 
behaves similarly to the current (Fig. pT] inset). For the 
case of only n = and n = 1 electron states the variance 
of the electron occupation number reads 

Var(n) = (n') - (n)' = (n) (1 - (n)) . (37) 

The increase of (n) > 1/2 signals a suppression of the 
fluctuations of the electronic level population analo- 
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FIG. 11: (Color online) (a) Current / (units eFo) as a function 
of V (units cjo/e) for different values of 7. (b) Ratio between 
the current obtained with the GME (/) and the current in the 
RWA approximation (/rwa) as a function of 7, at eV = IScJo- 
Inset: average occupation (n) of the electronic level. Other 



parameters: rig = 1/2, r = 9, A = 2 and ujc 



10^ 



gously to the fluctuations of the vibronic part. 



range, extends from the very fast vibrations ujq Fq 
to the slow, coherent regime where the off-diagonal 
elements of the reduced density matrix between energy 
eigenstates are paramount. 

In the coherent regime, two peculiar behaviors have 
been found. For intermediate frequencies, uq ^ Fq, 
the system can be described in terms of an effective 
thermal distribution with a temperature lower than 
that of the environment. The cooling phenomenon is 
accompanied by a decrease of position and momentum 
quadratures. For still slower oscillations, a strongly 
coherent regime is entered characterized by non classical 
behavior. A benchmark of this regime is a marked 
suppression of the vibron Fano factor, which can even 
attain suh-Poissonian values. 

This work is one of the first steps towards the understand- 
ing of dynamical properties of nano-electromechanical 
systems in the coherent regime and represents a rather 
tough numerical challenge, due to the slow convergence 
of the master equation solution and the need of a 
large number of basis states. Future investigations are 
certainly in order, to explore the regime of very strong 
electron-vibron coupling and the crossover towards the 
semi-classical regime. Further interesting issues to be 
investigated are higher order electronic properties, like 
the current fluctuations and their connection with the 
fluctuations of the mechanical part. 



IV. CONCLUSIONS 

In the present article we derived a generalized master 
equation to explore the steady-state properties of a 
nano-electromechanical system in a wide parameters 
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